
!------------------------------------------------------------------------!
!  The Community Multiscale Air Quality (CMAQ) system software is in     !
!  continuous development by various groups and is based on information  !
!  from these groups: Federal Government employees, contractors working  !
!  within a United States Government contract, and non-Federal sources   !
!  including research institutions.  These groups give the Government    !
!  permission to use, prepare derivative works of, and distribute copies !
!  of their work in the CMAQ system to the public and to permit others   !
!  to do so.  The United States Environmental Protection Agency          !
!  therefore grants similar permission to use the CMAQ system software,  !
!  but users are requested to provide copies of derivative works or      !
!  products designed to operate in the CMAQ system to the United States  !
!  Government without restrictions as to use by others.  Software        !
!  that is used with the CMAQ system but distributed under the GNU       !
!  General Public License or the GNU Lesser General Public License is    !
!  subject to their copyright restrictions.                              !
!------------------------------------------------------------------------!
 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! Linear Algebra Data and Routines File
! 
! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
!       (http://www.cs.vt.edu/~asandu/Software/KPP)
! KPP is distributed under GPL, the general public licence
!       (http://www.gnu.org/copyleft/gpl.html)
! (C) 1995-1997, V. Damian & A. Sandu, CGRER, Univ. Iowa
! (C) 1997-2005, A. Sandu, Michigan Tech, Virginia Tech
!     With important contributions from:
!        M. Damian, Villanova University, USA
!        R. Sander, Max-Planck Institute for Chemistry, Mainz, Germany
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



MODULE aqchem_LinearAlgebra

  USE aqchem_Parameters
  USE aqchem_JacobianSP

  IMPLICIT NONE

CONTAINS


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! SPARSE_UTIL - SPARSE utility functions
!   Arguments :
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppDecomp( JVS, IER )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!        Sparse LU factorization
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE aqchem_Parameters
  USE aqchem_JacobianSP

      INTEGER  :: IER
      REAL(kind=dp) :: JVS(LU_NONZERO), W(NVAR), a
      INTEGER  :: k, kk, j, jj

      a = 0. ! mz_rs_20050606
      IER = 0
      DO k=1,NVAR
        ! mz_rs_20050606: don't check if real value == 0
        ! IF ( JVS( LU_DIAG(k) ) .EQ. 0. ) THEN
        IF ( ABS(JVS(LU_DIAG(k))) < TINY(a) ) THEN
            IER = k
            RETURN
        END IF
        DO kk = LU_CROW(k), LU_CROW(k+1)-1
              W( LU_ICOL(kk) ) = JVS(kk)
        END DO
        DO kk = LU_CROW(k), LU_DIAG(k)-1
            j = LU_ICOL(kk)
            a = -W(j) / JVS( LU_DIAG(j) )
            W(j) = -a
            DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
               W( LU_ICOL(jj) ) = W( LU_ICOL(jj) ) + a*JVS(jj)
            END DO
         END DO
         DO kk = LU_CROW(k), LU_CROW(k+1)-1
            JVS(kk) = W( LU_ICOL(kk) )
         END DO
      END DO
      
END SUBROUTINE KppDecomp


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppDecompCmplx( JVS, IER )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!        Sparse LU factorization, complex
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE aqchem_Parameters
  USE aqchem_JacobianSP

      INTEGER        :: IER
      DOUBLE COMPLEX :: JVS(LU_NONZERO), W(NVAR), a
      REAL(kind=dp)  :: b = 0.0
      INTEGER        :: k, kk, j, jj

      IER = 0
      DO k=1,NVAR
        IF ( ABS(JVS(LU_DIAG(k))) < TINY(b) ) THEN
            IER = k
            RETURN
        END IF
        DO kk = LU_CROW(k), LU_CROW(k+1)-1
              W( LU_ICOL(kk) ) = JVS(kk)
        END DO
        DO kk = LU_CROW(k), LU_DIAG(k)-1
            j = LU_ICOL(kk)
            a = -W(j) / JVS( LU_DIAG(j) )
            W(j) = -a
            DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
               W( LU_ICOL(jj) ) = W( LU_ICOL(jj) ) + a*JVS(jj)
            END DO
         END DO
         DO kk = LU_CROW(k), LU_CROW(k+1)-1
            JVS(kk) = W( LU_ICOL(kk) )
         END DO
      END DO
      
END SUBROUTINE KppDecompCmplx


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppDecompCmplxR( JVSR, JVSI, IER )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Sparse LU factorization, complex
!   (Real and Imaginary parts are used instead of complex data type)     
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE aqchem_Parameters
  USE aqchem_JacobianSP

      INTEGER       :: IER
      REAL(kind=dp) :: JVSR(LU_NONZERO), JVSI(LU_NONZERO) 
      REAL(kind=dp) :: WR(NVAR), WI(NVAR), ar, ai, den
      INTEGER       :: k, kk, j, jj

      IER = 0
      ar  = 0.0
      DO k=1,NVAR
        IF (  ( ABS(JVSR(LU_DIAG(k))) < TINY(ar) ) .AND. &
              ( ABS(JVSI(LU_DIAG(k))) < TINY(ar) ) )  THEN
            IER = k
            RETURN
        END IF
        DO kk = LU_CROW(k), LU_CROW(k+1)-1
              WR( LU_ICOL(kk) ) = JVSR(kk)
              WI( LU_ICOL(kk) ) = JVSI(kk)
        END DO
        DO kk = LU_CROW(k), LU_DIAG(k)-1
            j = LU_ICOL(kk)
            den = JVSR(LU_DIAG(j))**2 + JVSI(LU_DIAG(j))**2
            ar = -(WR(j)*JVSR(LU_DIAG(j)) + WI(j)*JVSI(LU_DIAG(j)))/den
            ai = -(WI(j)*JVSR(LU_DIAG(j)) - WR(j)*JVSI(LU_DIAG(j)))/den
            WR(j) = -ar
            WI(j) = -ai
            DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
               WR( LU_ICOL(jj) ) = WR( LU_ICOL(jj) ) + ar*JVSR(jj) - ai*JVSI(jj)
               WI( LU_ICOL(jj) ) = WI( LU_ICOL(jj) ) + ar*JVSI(jj) + ai*JVSR(jj)
            END DO
         END DO
         DO kk = LU_CROW(k), LU_CROW(k+1)-1
            JVSR(kk) = WR( LU_ICOL(kk) )
            JVSI(kk) = WI( LU_ICOL(kk) )
         END DO
      END DO

END SUBROUTINE KppDecompCmplxR


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppSolveIndirect( JVS, X )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!        Sparse solve subroutine using indirect addressing
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE aqchem_Parameters
  USE aqchem_JacobianSP

      INTEGER  :: i, j
      REAL(kind=dp) :: JVS(LU_NONZERO), X(NVAR), sum

      DO i=1,NVAR
         DO j = LU_CROW(i), LU_DIAG(i)-1 
             X(i) = X(i) - JVS(j)*X(LU_ICOL(j));
         END DO  
      END DO

      DO i=NVAR,1,-1
        sum = X(i);
        DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
          sum = sum - JVS(j)*X(LU_ICOL(j));
        END DO
        X(i) = sum/JVS(LU_DIAG(i));
      END DO
      
END SUBROUTINE KppSolveIndirect


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppSolveTRIndirect( JVS, X )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!        Complex sparse solve transpose subroutine using indirect addressing
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE aqchem_Parameters
  USE aqchem_JacobianSP

      INTEGER       :: i, j
      REAL(kind=dp) :: JVS(LU_NONZERO), X(NVAR)

      DO i=1,NVAR
        X(i) = X(i)/JVS(LU_DIAG(i))
        ! subtract all nonzero elements in row i of JVS from X
        DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1
          X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i)
        END DO
      END DO

      DO i=NVAR, 1, -1
        ! subtract all nonzero elements in row i of JVS from X
        DO j=LU_CROW(i),LU_DIAG(i)-1
          X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i)
        END DO
      END DO
      
END SUBROUTINE KppSolveTRIndirect


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppSolveCmplx( JVS, X )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!        Complex sparse solve subroutine using indirect addressing
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE aqchem_Parameters
  USE aqchem_JacobianSP

      INTEGER        :: i, j
      DOUBLE COMPLEX :: JVS(LU_NONZERO), X(NVAR), sum

      DO i=1,NVAR
         DO j = LU_CROW(i), LU_DIAG(i)-1 
             X(i) = X(i) - JVS(j)*X(LU_ICOL(j));
         END DO  
      END DO

      DO i=NVAR,1,-1
        sum = X(i);
        DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
          sum = sum - JVS(j)*X(LU_ICOL(j));
        END DO
        X(i) = sum/JVS(LU_DIAG(i));
      END DO
      
END SUBROUTINE KppSolveCmplx

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppSolveCmplxR( JVSR, JVSI, XR, XI )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!   Complex sparse solve subroutine using indirect addressing
!   (Real and Imaginary parts are used instead of complex data type)     
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE aqchem_Parameters
  USE aqchem_JacobianSP

      INTEGER       ::  i, j
      REAL(kind=dp) ::  JVSR(LU_NONZERO), JVSI(LU_NONZERO), XR(NVAR), XI(NVAR), sumr, sumi, den

      DO i=1,NVAR
         DO j = LU_CROW(i), LU_DIAG(i)-1 
             XR(i) = XR(i) - (JVSR(j)*XR(LU_ICOL(j)) - JVSI(j)*XI(LU_ICOL(j)))
             XI(i) = XI(i) - (JVSR(j)*XI(LU_ICOL(j)) + JVSI(j)*XR(LU_ICOL(j)))
         END DO  
      END DO

      DO i=NVAR,1,-1
        sumr = XR(i); sumi = XI(i)
        DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
            sumr = sumr - (JVSR(j)*XR(LU_ICOL(j)) - JVSI(j)*XI(LU_ICOL(j)))
            sumi = sumi - (JVSR(j)*XI(LU_ICOL(j)) + JVSI(j)*XR(LU_ICOL(j)))
        END DO
        den   = JVSR(LU_DIAG(i))**2 + JVSI(LU_DIAG(i))**2
        XR(i) = (sumr*JVSR(LU_DIAG(i)) + sumi*JVSI(LU_DIAG(i)))/den
        XI(i) = (sumi*JVSR(LU_DIAG(i)) - sumr*JVSI(LU_DIAG(i)))/den
      END DO
      
END SUBROUTINE KppSolveCmplxR


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppSolveTRCmplx( JVS, X )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!        Complex sparse solve transpose subroutine using indirect addressing
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE aqchem_Parameters
  USE aqchem_JacobianSP

      INTEGER        :: i, j
      DOUBLE COMPLEX :: JVS(LU_NONZERO), X(NVAR)

      DO i=1,NVAR
        X(i) = X(i)/JVS(LU_DIAG(i))
        ! subtract all nonzero elements in row i of JVS from X
        DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1
          X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i)
        END DO
      END DO

      DO i=NVAR, 1, -1
        ! subtract all nonzero elements in row i of JVS from X
        DO j=LU_CROW(i),LU_DIAG(i)-1
          X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i)
        END DO
      END DO
      
END SUBROUTINE KppSolveTRCmplx


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE KppSolveTRCmplxR( JVSR, JVSI, XR, XI )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!   Complex sparse solve transpose subroutine using indirect addressing
!   (Real and Imaginary parts are used instead of complex data type)     
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE aqchem_Parameters
  USE aqchem_JacobianSP

      INTEGER       ::  i, j
      REAL(kind=dp) ::  JVSR(LU_NONZERO), JVSI(LU_NONZERO), XR(NVAR), XI(NVAR), den

      DO i=1,NVAR
        den   = JVSR(LU_DIAG(i))**2 + JVSI(LU_DIAG(i))**2
        XR(i) = (XR(i)*JVSR(LU_DIAG(i)) + XI(i)*JVSI(LU_DIAG(i)))/den
        XI(i) = (XI(i)*JVSR(LU_DIAG(i)) - XR(i)*JVSI(LU_DIAG(i)))/den
        ! subtract all nonzero elements in row i of JVS from X
        DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1
          XR(LU_ICOL(j)) = XR(LU_ICOL(j))-(JVSR(j)*XR(i) - JVSI(j)*XI(i))
          XI(LU_ICOL(j)) = XI(LU_ICOL(j))-(JVSI(j)*XR(i) + JVSR(j)*XI(i))
        END DO
      END DO

      DO i=NVAR, 1, -1
        ! subtract all nonzero elements in row i of JVS from X
        DO j=LU_CROW(i),LU_DIAG(i)-1
          XR(LU_ICOL(j)) = XR(LU_ICOL(j))-(JVSR(j)*XR(i) - JVSI(j)*XI(i))
          XI(LU_ICOL(j)) = XI(LU_ICOL(j))-(JVSI(j)*XR(i) + JVSR(j)*XI(i))
        END DO
      END DO
      
END SUBROUTINE KppSolveTRCmplxR


!
! Next few commented subroutines perform sparse big linear algebra
!
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!SUBROUTINE KppDecompBig( JVS, IP, IER )
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!!        Sparse LU factorization
!!        for the Runge Kutta (3n)x(3n) linear system
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!  USE aqchem_Parameters
!  USE aqchem_JacobianSP
!
!      INTEGER  :: IP3(3), IER, IP(3,NVAR)
!      REAL(kind=dp) :: JVS(3,3,LU_NONZERO), W(3,3,NVAR), a(3,3), E(3,3)
!      INTEGER  :: k, kk, j, jj
!
!      a = 0.0d0
!      IER = 0
!      DO k=1,NVAR
!        DO kk = LU_CROW(k), LU_CROW(k+1)-1
!              W( 1:3,1:3,LU_ICOL(kk) ) = JVS(1:3,1:3,kk)
!        END DO
!        DO kk = LU_CROW(k), LU_DIAG(k)-1
!            j = LU_ICOL(kk)
!            E(1:3,1:3) = JVS( 1:3,1:3,LU_DIAG(j) )
!            ! CALL DGETRF(3,3,E,3,IP3,IER) 
!            CALL FAC3(E,IP3,IER)
!            IF ( IER /= 0 )  RETURN
!            ! a = W(j) / JVS( LU_DIAG(j) )
!            a(1:3,1:3) = W( 1:3,1:3,j )
!            ! CALL DGETRS ('N',3,3,E,3,IP3,a,3,IER) 
!            CALL SOL3('N',E,IP3,a(1,1))
!            CALL SOL3('N',E,IP3,a(1,2))
!            CALL SOL3('N',E,IP3,a(1,3))
!            W(1:3,1:3,j) = a(1:3,1:3)
!            DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
!               W( 1:3,1:3,LU_ICOL(jj) ) = W( 1:3,1:3,LU_ICOL(jj) ) &
!                        - MATMUL( a(1:3,1:3) , JVS(1:3,1:3,jj) )
!            END DO
!         END DO
!         DO kk = LU_CROW(k), LU_CROW(k+1)-1
!            JVS(1:3,1:3,kk) = W( 1:3,1:3,LU_ICOL(kk) )
!         END DO
!      END DO
!
!      DO k=1,NVAR
!         ! CALL WGEFA(JVS(1,1,LU_DIAG(k)),3,3,IP(1,k),IER)
!         ! CALL DGETRF(3,3,JVS(1,1,LU_DIAG(k)),3,IP(1,k),IER)
!         CALL FAC3(JVS(1,1,LU_DIAG(k)),IP(1,k),IER)
!         IF ( IER /= 0 )  RETURN
!      END DO 
!      
!END SUBROUTINE KppDecompBig
!
!
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!SUBROUTINE KppSolveBig( JVS, IP, X )
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!!        Sparse solve subroutine using indirect addressing
!!        for the Runge Kutta (3n)x(3n) linear system
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!  USE aqchem_Parameters
!  USE aqchem_JacobianSP
!
!      INTEGER  :: i, j, k, m, IP3(3), IP(3,NVAR), IER
!      REAL(kind=dp) :: JVS(3,3,LU_NONZERO), X(3,NVAR), sum(3)
!
!      DO i=1,NVAR
!        DO j = LU_CROW(i), LU_DIAG(i)-1 
!          !X(1:3,i) = X(1:3,i) - MATMUL(JVS(1:3,1:3,j),X(1:3,LU_ICOL(j)));
!          DO k=1,3
!            DO m=1,3
!               X(k,i) = X(k,i) - JVS(k,m,j)*X(m,LU_ICOL(j))
!            END DO
!          END DO
!        END DO  
!      END DO
!
!      DO i=NVAR,1,-1
!        sum(1:3) = X(1:3,i);
!        DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
!          !sum(1:3) = sum(1:3) - MATMUL(JVS(1:3,1:3,j),X(1:3,LU_ICOL(j)));
!          DO k=1,3
!            DO m=1,3
!               sum(k) = sum(k) - JVS(k,m,j)*X(m,LU_ICOL(j))
!            END DO
!          END DO
!        END DO
!        ! X(i) = sum/JVS(LU_DIAG(i));
!        ! CALL DGETRS ('N',3,1,JVS(1:3,1:3,LU_DIAG(i)),3,IP(1,i),sum,3,0) 
!        ! CALL WGESL('N',JVS(1,1,LU_DIAG(i)),3,3,IP(1,i),sum)
!        CALL SOL3('N',JVS(1,1,LU_DIAG(i)),IP(1,i),sum)
!        X(1:3,i) = sum(1:3)
!      END DO
!      
!END SUBROUTINE KppSolveBig
!
!
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!SUBROUTINE KppSolveBigTR( JVS, IP, X )
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!!        Big sparse transpose solve using indirect addressing
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!  USE aqchem_Parameters
!  USE aqchem_JacobianSP
!
!      INTEGER       :: i, j, k, m, IP(3,NVAR)
!      REAL(kind=dp) :: JVS(3,3,LU_NONZERO), X(3,NVAR)
!
!      DO i=1,NVAR
!        ! X(i) = X(i)/JVS(LU_DIAG(i))
!        CALL SOL3('T',JVS(1,1,LU_DIAG(i)),IP(1,i),X(1,i))
!        DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1
!          !X(1:3,LU_ICOL(j)) = X(1:3,LU_ICOL(j)) &
!          !    - MATMUL( TRANSPOSE(JVS(1:3,1:3,j)), X(1:3,i) )
!          DO k=1,3
!            DO m=1,3
!               X(k,LU_ICOL(j)) = X(k,LU_ICOL(j)) - JVS(m,k,j)*X(m,i)
!            END DO
!          END DO
!        END DO
!      END DO
!
!      DO i=NVAR, 1, -1
!        DO j=LU_CROW(i),LU_DIAG(i)-1
!          !X(1:3,LU_ICOL(j)) = X(1:3,LU_ICOL(j)) &
!          !   - MATMUL( TRANSPOSE(JVS(1:3,1:3,j)), X(1:3,i) )
!          DO k=1,3
!            DO m=1,3
!               X(k,LU_ICOL(j)) = X(k,LU_ICOL(j)) - JVS(m,k,j)*X(m,i)
!            END DO
!          END DO
!        END DO
!      END DO
!      
!END SUBROUTINE KppSolveBigTR
!
!
!
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!SUBROUTINE FAC3(A,IPVT,INFO)
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!!     FAC3 FACTORS THE MATRIX A (3,3) BY
!!           GAUSS ELIMINATION WITH PARTIAL PIVOTING
!!     LINPACK - LIKE 
!!
!!     Remove comments to perform pivoting
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!!
!      REAL(kind=dp) :: A(3,3)
!      INTEGER       :: IPVT(3),INFO
!!      INTEGER       :: L
!!      REAL(kind=dp) :: t, dmax, da, TMP(3)
!      REAL(kind=dp), PARAMETER :: ZERO = 0.0, ONE = 1.0
!
!      info = 0
!!      t = TINY(da)
!!      
!!      da = ABS(A(1,1)); L = 1
!!      IF ( ABS(A(2,1))>da ) THEN
!!        da = ABS(A(2,1)); L = 2
!!        IF ( ABS(A(3,1))>da ) THEN
!!          L = 3
!!        END IF  
!!      END IF  
!!      IPVT(1)  = L
!!      IF (L /=1 ) THEN
!!         TMP(1:3) = A(L,1:3)
!!         A(L,1:3) = A(1,1:3)
!!         A(1,1:3) = TMP(1:3)
!!      END IF
!!      IF (ABS(A(1,1)) < t) THEN
!!         info = 1
!!         return
!!      END IF   
!!
!      A(2,1) = A(2,1)/A(1,1)
!      A(2,2) = A(2,2) - A(2,1)*A(1,2)
!      A(2,3) = A(2,3) - A(2,1)*A(1,3)
!      A(3,1) = A(3,1)/A(1,1)
!      A(3,2) = A(3,2) - A(3,1)*A(1,2)
!      A(3,3) = A(3,3) - A(3,1)*A(1,3)
!      
!!      IPVT(2)  = 2
!!      IF (ABS(A(3,2))>ABS(A(2,2))) THEN
!!         IPVT(2)  = 3
!!         TMP(2:3) = A(3,2:3)
!!         A(3,2:3) = A(2,2:3)
!!         A(2,2:3) = TMP(2:3)
!!      END IF
!!      IF (ABS(A(2,2)) < t) THEN
!!         info = 1
!!         return
!!      END IF   
!!      
!      A(3,2)   = A(3,2)/A(2,2)
!      A(3,3)   = A(3,3) - A(3,2)*A(2,3)
!      IPVT(3)  = 3
!      
!END SUBROUTINE FAC3
!
!
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!SUBROUTINE SOL3(Trans,A,IPVT,b)
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!!     SOL3 solves the system 3x3
!!     A * x = b  or  trans(a) * x = b
!!     using the factors computed by WGEFA.
!!
!!     Trans      = 'N'   to solve  A*x = b ,
!!                = 'T'   to solve  transpose(A)*x = b
!!     LINPACK - LIKE 
!!
!!     Remove comments to use pivoting
!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!      CHARACTER     :: Trans
!      REAL(kind=dp) :: a(3,3),b(3)
!      INTEGER       :: IPVT(3)
!!      INTEGER       :: L
!!      REAL(kind=dp) :: TMP
!      
!      SELECT CASE (Trans)
!
!      CASE ('n','N')  !  Solve  A * x = b
!
!!     Solve  L*y = b
!!         L = IPVT(1)
!!         IF (L /= 1) THEN
!!            TMP = B(1); B(1) = B(L); B(L) = TMP
!!         END IF
!         b(2) = b(2)-A(2,1)*b(1)
!         b(3) = b(3)-A(3,1)*b(1)
!         
!!         L = IPVT(2)
!!         IF (L /= 2) THEN
!!            TMP = B(2); B(2) = B(L); B(L) = TMP
!!         END IF
!         b(3) = b(3)-A(3,2)*b(2)
!
!!     Solve  U*x = y
!         b(3) = b(3)/A(3,3)
!         b(2) = (b(2)-A(2,3)*b(3))/A(2,2)
!         b(1) = (b(1)-A(1,3)*b(3)-A(1,2)*b(2))/A(1,1)
!      
!      
!      CASE ('t','T')  !  Solve transpose(A) * x = b
!
!!      Solve transpose(U)*y = b
!         b(1) = b(1)/A(1,1)
!         b(2) = (b(2)-A(1,2)*b(1))/A(2,2)
!         b(3) = (b(3)-A(1,3)*b(1)-A(2,3)*b(2))/A(3,3)
!
!!      Solve transpose(L)*x = y
!         b(2) = b(2)-A(3,2)*b(3)
!!         L = ipvt(2)
!!         IF (L /= 2) THEN
!!            TMP = B(2); B(2) = B(L); B(L) = TMP
!!         END IF
!         b(1) = b(1)-A(3,1)*b(3)-A(2,1)*b(2)
!!         L = ipvt(1)
!!         IF (L /= 1) THEN
!!            TMP = B(1); B(1) = B(L); B(L) = TMP
!!         END IF
!   
!      END SELECT
!
!END SUBROUTINE SOL3

! End of SPARSE_UTIL function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! KppSolve - sparse back substitution
!   Arguments :
!      JVS       - sparse Jacobian of variables
!      X         - Vector for variables
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SUBROUTINE KppSolve ( JVS, X )

! JVS - sparse Jacobian of variables
  REAL(kind=dp) :: JVS(LU_NONZERO)
! X - Vector for variables
  REAL(kind=dp) :: X(NVAR)

  X(9) = X(9)-JVS(9)*X(8)
  X(10) = X(10)-JVS(11)*X(6)
  X(11) = X(11)-JVS(13)*X(7)
  X(17) = X(17)-JVS(20)*X(4)
  X(34) = X(34)-JVS(61)*X(9)
  X(35) = X(35)-JVS(63)*X(12)
  X(36) = X(36)-JVS(65)*X(13)
  X(37) = X(37)-JVS(67)*X(10)
  X(40) = X(40)-JVS(74)*X(11)
  X(42) = X(42)-JVS(78)*X(14)
  X(43) = X(43)-JVS(80)*X(15)
  X(44) = X(44)-JVS(82)*X(16)
  X(45) = X(45)-JVS(84)*X(17)
  X(54) = X(54)-JVS(102)*X(53)
  X(56) = X(56)-JVS(106)*X(55)
  X(60) = X(60)-JVS(114)*X(46)
  X(61) = X(61)-JVS(118)*X(1)-JVS(119)*X(60)
  X(62) = X(62)-JVS(122)*X(48)
  X(63) = X(63)-JVS(126)*X(2)-JVS(127)*X(62)
  X(64) = X(64)-JVS(130)*X(62)-JVS(131)*X(63)
  X(65) = X(65)-JVS(134)*X(51)
  X(66) = X(66)-JVS(138)*X(65)
  X(68) = X(68)-JVS(144)*X(58)
  X(69) = X(69)-JVS(148)*X(3)-JVS(149)*X(68)
  X(71) = X(71)-JVS(155)*X(47)
  X(72) = X(72)-JVS(159)*X(70)-JVS(160)*X(71)
  X(73) = X(73)-JVS(163)*X(67)
  X(74) = X(74)-JVS(167)*X(49)
  X(75) = X(75)-JVS(171)*X(57)
  X(76) = X(76)-JVS(175)*X(5)-JVS(176)*X(12)-JVS(177)*X(13)-JVS(178)*X(73)-JVS(179)*X(74)-JVS(180)*X(75)
  X(77) = X(77)-JVS(188)*X(52)
  X(78) = X(78)-JVS(192)*X(12)-JVS(193)*X(13)
  X(79) = X(79)-JVS(198)*X(12)-JVS(199)*X(13)-JVS(200)*X(60)-JVS(201)*X(61)-JVS(202)*X(64)-JVS(203)*X(65)-JVS(204)*X(66)&
            &-JVS(205)*X(67)-JVS(206)*X(68)-JVS(207)*X(69)-JVS(208)*X(70)-JVS(209)*X(71)-JVS(210)*X(72)-JVS(211)*X(73)&
            &-JVS(212)*X(74)-JVS(213)*X(75)-JVS(214)*X(76)-JVS(215)*X(77)-JVS(216)*X(78)
  X(80) = X(80)-JVS(221)*X(50)-JVS(222)*X(78)-JVS(223)*X(79)
  X(81) = X(81)-JVS(227)*X(12)-JVS(228)*X(13)-JVS(229)*X(74)-JVS(230)*X(75)-JVS(231)*X(77)-JVS(232)*X(78)-JVS(233)*X(79)&
            &-JVS(234)*X(80)
  X(82) = X(82)-JVS(237)*X(12)-JVS(238)*X(13)-JVS(239)*X(59)-JVS(240)*X(79)-JVS(241)*X(80)-JVS(242)*X(81)
  X(82) = X(82)/JVS(243)
  X(81) = (X(81)-JVS(236)*X(82))/(JVS(235))
  X(80) = (X(80)-JVS(225)*X(81)-JVS(226)*X(82))/(JVS(224))
  X(79) = (X(79)-JVS(218)*X(80)-JVS(219)*X(81)-JVS(220)*X(82))/(JVS(217))
  X(78) = (X(78)-JVS(195)*X(79)-JVS(196)*X(80)-JVS(197)*X(81))/(JVS(194))
  X(77) = (X(77)-JVS(190)*X(79)-JVS(191)*X(81))/(JVS(189))
  X(76) = (X(76)-JVS(182)*X(77)-JVS(183)*X(78)-JVS(184)*X(79)-JVS(185)*X(80)-JVS(186)*X(81)-JVS(187)*X(82))/(JVS(181))
  X(75) = (X(75)-JVS(173)*X(79)-JVS(174)*X(81))/(JVS(172))
  X(74) = (X(74)-JVS(169)*X(79)-JVS(170)*X(81))/(JVS(168))
  X(73) = (X(73)-JVS(165)*X(76)-JVS(166)*X(79))/(JVS(164))
  X(72) = (X(72)-JVS(162)*X(79))/(JVS(161))
  X(71) = (X(71)-JVS(157)*X(72)-JVS(158)*X(79))/(JVS(156))
  X(70) = (X(70)-JVS(153)*X(72)-JVS(154)*X(79))/(JVS(152))
  X(69) = (X(69)-JVS(151)*X(79))/(JVS(150))
  X(68) = (X(68)-JVS(146)*X(69)-JVS(147)*X(79))/(JVS(145))
  X(67) = (X(67)-JVS(142)*X(73)-JVS(143)*X(79))/(JVS(141))
  X(66) = (X(66)-JVS(140)*X(79))/(JVS(139))
  X(65) = (X(65)-JVS(136)*X(66)-JVS(137)*X(79))/(JVS(135))
  X(64) = (X(64)-JVS(133)*X(79))/(JVS(132))
  X(63) = (X(63)-JVS(129)*X(64))/(JVS(128))
  X(62) = (X(62)-JVS(124)*X(63)-JVS(125)*X(64))/(JVS(123))
  X(61) = (X(61)-JVS(121)*X(79))/(JVS(120))
  X(60) = (X(60)-JVS(116)*X(61)-JVS(117)*X(79))/(JVS(115))
  X(59) = (X(59)-JVS(113)*X(82))/(JVS(112))
  X(58) = (X(58)-JVS(111)*X(68))/(JVS(110))
  X(57) = (X(57)-JVS(109)*X(75))/(JVS(108))
  X(56) = X(56)/JVS(107)
  X(55) = (X(55)-JVS(105)*X(56))/(JVS(104))
  X(54) = X(54)/JVS(103)
  X(53) = (X(53)-JVS(101)*X(54))/(JVS(100))
  X(52) = (X(52)-JVS(99)*X(77))/(JVS(98))
  X(51) = (X(51)-JVS(97)*X(65))/(JVS(96))
  X(50) = (X(50)-JVS(95)*X(80))/(JVS(94))
  X(49) = (X(49)-JVS(93)*X(74))/(JVS(92))
  X(48) = (X(48)-JVS(91)*X(62))/(JVS(90))
  X(47) = (X(47)-JVS(89)*X(71))/(JVS(88))
  X(46) = (X(46)-JVS(87)*X(60))/(JVS(86))
  X(45) = X(45)/JVS(85)
  X(44) = X(44)/JVS(83)
  X(43) = X(43)/JVS(81)
  X(42) = X(42)/JVS(79)
  X(41) = (X(41)-JVS(77)*X(79))/(JVS(76))
  X(40) = X(40)/JVS(75)
  X(39) = (X(39)-JVS(72)*X(53)-JVS(73)*X(55))/(JVS(71))
  X(38) = (X(38)-JVS(70)*X(39))/(JVS(69))
  X(37) = X(37)/JVS(68)
  X(36) = X(36)/JVS(66)
  X(35) = X(35)/JVS(64)
  X(34) = X(34)/JVS(62)
  X(33) = (X(33)-JVS(60)*X(69))/(JVS(59))
  X(32) = (X(32)-JVS(58)*X(63))/(JVS(57))
  X(31) = (X(31)-JVS(56)*X(61))/(JVS(55))
  X(30) = (X(30)-JVS(54)*X(55))/(JVS(53))
  X(29) = (X(29)-JVS(52)*X(53))/(JVS(51))
  X(28) = (X(28)-JVS(50)*X(68))/(JVS(49))
  X(27) = (X(27)-JVS(46)*X(67)-JVS(47)*X(73)-JVS(48)*X(76))/(JVS(45))
  X(26) = (X(26)-JVS(44)*X(75))/(JVS(43))
  X(25) = (X(25)-JVS(42)*X(77))/(JVS(41))
  X(24) = (X(24)-JVS(39)*X(65)-JVS(40)*X(66))/(JVS(38))
  X(23) = (X(23)-JVS(37)*X(80))/(JVS(36))
  X(22) = (X(22)-JVS(35)*X(74))/(JVS(34))
  X(21) = (X(21)-JVS(33)*X(62))/(JVS(32))
  X(20) = (X(20)-JVS(29)*X(70)-JVS(30)*X(71)-JVS(31)*X(72))/(JVS(28))
  X(19) = (X(19)-JVS(27)*X(60))/(JVS(26))
  X(18) = (X(18)-JVS(23)*X(78)-JVS(24)*X(81)-JVS(25)*X(82))/(JVS(22))
  X(17) = X(17)/JVS(21)
  X(16) = X(16)/JVS(19)
  X(15) = X(15)/JVS(18)
  X(14) = X(14)/JVS(17)
  X(13) = X(13)/JVS(16)
  X(12) = X(12)/JVS(15)
  X(11) = X(11)/JVS(14)
  X(10) = X(10)/JVS(12)
  X(9) = X(9)/JVS(10)
  X(8) = X(8)/JVS(8)
  X(7) = X(7)/JVS(7)
  X(6) = X(6)/JVS(6)
  X(5) = X(5)/JVS(5)
  X(4) = X(4)/JVS(4)
  X(3) = X(3)/JVS(3)
  X(2) = X(2)/JVS(2)
  X(1) = X(1)/JVS(1)
      
END SUBROUTINE KppSolve

! End of KppSolve function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! KppSolveTR - sparse, transposed back substitution
!   Arguments :
!      JVS       - sparse Jacobian of variables
!      X         - Vector for variables
!      XX        - Vector for output variables
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SUBROUTINE KppSolveTR ( JVS, X, XX )

! JVS - sparse Jacobian of variables
  REAL(kind=dp) :: JVS(LU_NONZERO)
! X - Vector for variables
  REAL(kind=dp) :: X(NVAR)
! XX - Vector for output variables
  REAL(kind=dp) :: XX(NVAR)

  XX(1) = X(1)/JVS(1)
  XX(2) = X(2)/JVS(2)
  XX(3) = X(3)/JVS(3)
  XX(4) = X(4)/JVS(4)
  XX(5) = X(5)/JVS(5)
  XX(6) = X(6)/JVS(6)
  XX(7) = X(7)/JVS(7)
  XX(8) = X(8)/JVS(8)
  XX(9) = X(9)/JVS(10)
  XX(10) = X(10)/JVS(12)
  XX(11) = X(11)/JVS(14)
  XX(12) = X(12)/JVS(15)
  XX(13) = X(13)/JVS(16)
  XX(14) = X(14)/JVS(17)
  XX(15) = X(15)/JVS(18)
  XX(16) = X(16)/JVS(19)
  XX(17) = X(17)/JVS(21)
  XX(18) = X(18)/JVS(22)
  XX(19) = X(19)/JVS(26)
  XX(20) = X(20)/JVS(28)
  XX(21) = X(21)/JVS(32)
  XX(22) = X(22)/JVS(34)
  XX(23) = X(23)/JVS(36)
  XX(24) = X(24)/JVS(38)
  XX(25) = X(25)/JVS(41)
  XX(26) = X(26)/JVS(43)
  XX(27) = X(27)/JVS(45)
  XX(28) = X(28)/JVS(49)
  XX(29) = X(29)/JVS(51)
  XX(30) = X(30)/JVS(53)
  XX(31) = X(31)/JVS(55)
  XX(32) = X(32)/JVS(57)
  XX(33) = X(33)/JVS(59)
  XX(34) = X(34)/JVS(62)
  XX(35) = X(35)/JVS(64)
  XX(36) = X(36)/JVS(66)
  XX(37) = X(37)/JVS(68)
  XX(38) = X(38)/JVS(69)
  XX(39) = (X(39)-JVS(70)*XX(38))/(JVS(71))
  XX(40) = X(40)/JVS(75)
  XX(41) = X(41)/JVS(76)
  XX(42) = X(42)/JVS(79)
  XX(43) = X(43)/JVS(81)
  XX(44) = X(44)/JVS(83)
  XX(45) = X(45)/JVS(85)
  XX(46) = X(46)/JVS(86)
  XX(47) = X(47)/JVS(88)
  XX(48) = X(48)/JVS(90)
  XX(49) = X(49)/JVS(92)
  XX(50) = X(50)/JVS(94)
  XX(51) = X(51)/JVS(96)
  XX(52) = X(52)/JVS(98)
  XX(53) = (X(53)-JVS(52)*XX(29)-JVS(72)*XX(39))/(JVS(100))
  XX(54) = (X(54)-JVS(101)*XX(53))/(JVS(103))
  XX(55) = (X(55)-JVS(54)*XX(30)-JVS(73)*XX(39))/(JVS(104))
  XX(56) = (X(56)-JVS(105)*XX(55))/(JVS(107))
  XX(57) = X(57)/JVS(108)
  XX(58) = X(58)/JVS(110)
  XX(59) = X(59)/JVS(112)
  XX(60) = (X(60)-JVS(27)*XX(19)-JVS(87)*XX(46))/(JVS(115))
  XX(61) = (X(61)-JVS(56)*XX(31)-JVS(116)*XX(60))/(JVS(120))
  XX(62) = (X(62)-JVS(33)*XX(21)-JVS(91)*XX(48))/(JVS(123))
  XX(63) = (X(63)-JVS(58)*XX(32)-JVS(124)*XX(62))/(JVS(128))
  XX(64) = (X(64)-JVS(125)*XX(62)-JVS(129)*XX(63))/(JVS(132))
  XX(65) = (X(65)-JVS(39)*XX(24)-JVS(97)*XX(51))/(JVS(135))
  XX(66) = (X(66)-JVS(40)*XX(24)-JVS(136)*XX(65))/(JVS(139))
  XX(67) = (X(67)-JVS(46)*XX(27))/(JVS(141))
  XX(68) = (X(68)-JVS(50)*XX(28)-JVS(111)*XX(58))/(JVS(145))
  XX(69) = (X(69)-JVS(60)*XX(33)-JVS(146)*XX(68))/(JVS(150))
  XX(70) = (X(70)-JVS(29)*XX(20))/(JVS(152))
  XX(71) = (X(71)-JVS(30)*XX(20)-JVS(89)*XX(47))/(JVS(156))
  XX(72) = (X(72)-JVS(31)*XX(20)-JVS(153)*XX(70)-JVS(157)*XX(71))/(JVS(161))
  XX(73) = (X(73)-JVS(47)*XX(27)-JVS(142)*XX(67))/(JVS(164))
  XX(74) = (X(74)-JVS(35)*XX(22)-JVS(93)*XX(49))/(JVS(168))
  XX(75) = (X(75)-JVS(44)*XX(26)-JVS(109)*XX(57))/(JVS(172))
  XX(76) = (X(76)-JVS(48)*XX(27)-JVS(165)*XX(73))/(JVS(181))
  XX(77) = (X(77)-JVS(42)*XX(25)-JVS(99)*XX(52)-JVS(182)*XX(76))/(JVS(189))
  XX(78) = (X(78)-JVS(23)*XX(18)-JVS(183)*XX(76))/(JVS(194))
  XX(79) = (X(79)-JVS(77)*XX(41)-JVS(117)*XX(60)-JVS(121)*XX(61)-JVS(133)*XX(64)-JVS(137)*XX(65)-JVS(140)*XX(66)&
             &-JVS(143)*XX(67)-JVS(147)*XX(68)-JVS(151)*XX(69)-JVS(154)*XX(70)-JVS(158)*XX(71)-JVS(162)*XX(72)-JVS(166)&
             &*XX(73)-JVS(169)*XX(74)-JVS(173)*XX(75)-JVS(184)*XX(76)-JVS(190)*XX(77)-JVS(195)*XX(78))/(JVS(217))
  XX(80) = (X(80)-JVS(37)*XX(23)-JVS(95)*XX(50)-JVS(185)*XX(76)-JVS(196)*XX(78)-JVS(218)*XX(79))/(JVS(224))
  XX(81) = (X(81)-JVS(24)*XX(18)-JVS(170)*XX(74)-JVS(174)*XX(75)-JVS(186)*XX(76)-JVS(191)*XX(77)-JVS(197)*XX(78)&
             &-JVS(219)*XX(79)-JVS(225)*XX(80))/(JVS(235))
  XX(82) = (X(82)-JVS(25)*XX(18)-JVS(113)*XX(59)-JVS(187)*XX(76)-JVS(220)*XX(79)-JVS(226)*XX(80)-JVS(236)*XX(81))&
             &/(JVS(243))
  XX(82) = XX(82)
  XX(81) = XX(81)-JVS(242)*XX(82)
  XX(80) = XX(80)-JVS(234)*XX(81)-JVS(241)*XX(82)
  XX(79) = XX(79)-JVS(223)*XX(80)-JVS(233)*XX(81)-JVS(240)*XX(82)
  XX(78) = XX(78)-JVS(216)*XX(79)-JVS(222)*XX(80)-JVS(232)*XX(81)
  XX(77) = XX(77)-JVS(215)*XX(79)-JVS(231)*XX(81)
  XX(76) = XX(76)-JVS(214)*XX(79)
  XX(75) = XX(75)-JVS(180)*XX(76)-JVS(213)*XX(79)-JVS(230)*XX(81)
  XX(74) = XX(74)-JVS(179)*XX(76)-JVS(212)*XX(79)-JVS(229)*XX(81)
  XX(73) = XX(73)-JVS(178)*XX(76)-JVS(211)*XX(79)
  XX(72) = XX(72)-JVS(210)*XX(79)
  XX(71) = XX(71)-JVS(160)*XX(72)-JVS(209)*XX(79)
  XX(70) = XX(70)-JVS(159)*XX(72)-JVS(208)*XX(79)
  XX(69) = XX(69)-JVS(207)*XX(79)
  XX(68) = XX(68)-JVS(149)*XX(69)-JVS(206)*XX(79)
  XX(67) = XX(67)-JVS(163)*XX(73)-JVS(205)*XX(79)
  XX(66) = XX(66)-JVS(204)*XX(79)
  XX(65) = XX(65)-JVS(138)*XX(66)-JVS(203)*XX(79)
  XX(64) = XX(64)-JVS(202)*XX(79)
  XX(63) = XX(63)-JVS(131)*XX(64)
  XX(62) = XX(62)-JVS(127)*XX(63)-JVS(130)*XX(64)
  XX(61) = XX(61)-JVS(201)*XX(79)
  XX(60) = XX(60)-JVS(119)*XX(61)-JVS(200)*XX(79)
  XX(59) = XX(59)-JVS(239)*XX(82)
  XX(58) = XX(58)-JVS(144)*XX(68)
  XX(57) = XX(57)-JVS(171)*XX(75)
  XX(56) = XX(56)
  XX(55) = XX(55)-JVS(106)*XX(56)
  XX(54) = XX(54)
  XX(53) = XX(53)-JVS(102)*XX(54)
  XX(52) = XX(52)-JVS(188)*XX(77)
  XX(51) = XX(51)-JVS(134)*XX(65)
  XX(50) = XX(50)-JVS(221)*XX(80)
  XX(49) = XX(49)-JVS(167)*XX(74)
  XX(48) = XX(48)-JVS(122)*XX(62)
  XX(47) = XX(47)-JVS(155)*XX(71)
  XX(46) = XX(46)-JVS(114)*XX(60)
  XX(45) = XX(45)
  XX(44) = XX(44)
  XX(43) = XX(43)
  XX(42) = XX(42)
  XX(41) = XX(41)
  XX(40) = XX(40)
  XX(39) = XX(39)
  XX(38) = XX(38)
  XX(37) = XX(37)
  XX(36) = XX(36)
  XX(35) = XX(35)
  XX(34) = XX(34)
  XX(33) = XX(33)
  XX(32) = XX(32)
  XX(31) = XX(31)
  XX(30) = XX(30)
  XX(29) = XX(29)
  XX(28) = XX(28)
  XX(27) = XX(27)
  XX(26) = XX(26)
  XX(25) = XX(25)
  XX(24) = XX(24)
  XX(23) = XX(23)
  XX(22) = XX(22)
  XX(21) = XX(21)
  XX(20) = XX(20)
  XX(19) = XX(19)
  XX(18) = XX(18)
  XX(17) = XX(17)-JVS(84)*XX(45)
  XX(16) = XX(16)-JVS(82)*XX(44)
  XX(15) = XX(15)-JVS(80)*XX(43)
  XX(14) = XX(14)-JVS(78)*XX(42)
  XX(13) = XX(13)-JVS(65)*XX(36)-JVS(177)*XX(76)-JVS(193)*XX(78)-JVS(199)*XX(79)-JVS(228)*XX(81)-JVS(238)*XX(82)
  XX(12) = XX(12)-JVS(63)*XX(35)-JVS(176)*XX(76)-JVS(192)*XX(78)-JVS(198)*XX(79)-JVS(227)*XX(81)-JVS(237)*XX(82)
  XX(11) = XX(11)-JVS(74)*XX(40)
  XX(10) = XX(10)-JVS(67)*XX(37)
  XX(9) = XX(9)-JVS(61)*XX(34)
  XX(8) = XX(8)-JVS(9)*XX(9)
  XX(7) = XX(7)-JVS(13)*XX(11)
  XX(6) = XX(6)-JVS(11)*XX(10)
  XX(5) = XX(5)-JVS(175)*XX(76)
  XX(4) = XX(4)-JVS(20)*XX(17)
  XX(3) = XX(3)-JVS(148)*XX(69)
  XX(2) = XX(2)-JVS(126)*XX(63)
  XX(1) = XX(1)-JVS(118)*XX(61)
      
END SUBROUTINE KppSolveTR

! End of KppSolveTR function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! BLAS_UTIL - BLAS-LIKE utility functions
!   Arguments :
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

!--------------------------------------------------------------
!
! BLAS/LAPACK-like subroutines used by the integration algorithms
! It is recommended to replace them by calls to the optimized
!      BLAS/LAPACK library for your machine
!
!  (C) Adrian Sandu, Aug. 2004
!      Virginia Polytechnic Institute and State University
!--------------------------------------------------------------


!--------------------------------------------------------------
      SUBROUTINE WCOPY(N,X,incX,Y,incY)
!--------------------------------------------------------------
!     copies a vector, x, to a vector, y:  y <- x
!     only for incX=incY=1
!     after BLAS
!     replace this by the function from the optimized BLAS implementation:
!         CALL  SCOPY(N,X,1,Y,1)   or   CALL  DCOPY(N,X,1,Y,1)
!--------------------------------------------------------------
!     USE aqchem_Precision
      
      INTEGER  :: i,incX,incY,M,MP1,N
      REAL(kind=dp) :: X(N),Y(N)

      IF (N.LE.0) RETURN

      M = MOD(N,8)
      IF( M .NE. 0 ) THEN
        DO i = 1,M
          Y(i) = X(i)
        END DO
        IF( N .LT. 8 ) RETURN
      END IF    
      MP1 = M+1
      DO i = MP1,N,8
        Y(i) = X(i)
        Y(i + 1) = X(i + 1)
        Y(i + 2) = X(i + 2)
        Y(i + 3) = X(i + 3)
        Y(i + 4) = X(i + 4)
        Y(i + 5) = X(i + 5)
        Y(i + 6) = X(i + 6)
        Y(i + 7) = X(i + 7)
      END DO

      END SUBROUTINE WCOPY


!--------------------------------------------------------------
      SUBROUTINE WAXPY(N,Alpha,X,incX,Y,incY)
!--------------------------------------------------------------
!     constant times a vector plus a vector: y <- y + Alpha*x
!     only for incX=incY=1
!     after BLAS
!     replace this by the function from the optimized BLAS implementation:
!         CALL SAXPY(N,Alpha,X,1,Y,1) or  CALL DAXPY(N,Alpha,X,1,Y,1)
!--------------------------------------------------------------

      INTEGER  :: i,incX,incY,M,MP1,N
      REAL(kind=dp) :: X(N),Y(N),Alpha
      REAL(kind=dp), PARAMETER :: ZERO = 0.0_dp

      IF (Alpha .EQ. ZERO) RETURN
      IF (N .LE. 0) RETURN

      M = MOD(N,4)
      IF( M .NE. 0 ) THEN
        DO i = 1,M
          Y(i) = Y(i) + Alpha*X(i)
        END DO
        IF( N .LT. 4 ) RETURN
      END IF
      MP1 = M + 1
      DO i = MP1,N,4
        Y(i) = Y(i) + Alpha*X(i)
        Y(i + 1) = Y(i + 1) + Alpha*X(i + 1)
        Y(i + 2) = Y(i + 2) + Alpha*X(i + 2)
        Y(i + 3) = Y(i + 3) + Alpha*X(i + 3)
      END DO
      
      END SUBROUTINE WAXPY



!--------------------------------------------------------------
      SUBROUTINE WSCAL(N,Alpha,X,incX)
!--------------------------------------------------------------
!     constant times a vector: x(1:N) <- Alpha*x(1:N) 
!     only for incX=incY=1
!     after BLAS
!     replace this by the function from the optimized BLAS implementation:
!         CALL SSCAL(N,Alpha,X,1) or  CALL DSCAL(N,Alpha,X,1)
!--------------------------------------------------------------

      INTEGER  :: i,incX,M,MP1,N
      REAL(kind=dp)  :: X(N),Alpha
      REAL(kind=dp), PARAMETER  :: ZERO=0.0_dp, ONE=1.0_dp

      IF (Alpha .EQ. ONE) RETURN
      IF (N .LE. 0) RETURN

      M = MOD(N,5)
      IF( M .NE. 0 ) THEN
        IF (Alpha .EQ. (-ONE)) THEN
          DO i = 1,M
            X(i) = -X(i)
          END DO
        ELSEIF (Alpha .EQ. ZERO) THEN
          DO i = 1,M
            X(i) = ZERO
          END DO
        ELSE
          DO i = 1,M
            X(i) = Alpha*X(i)
          END DO
        END IF
        IF( N .LT. 5 ) RETURN
      END IF
      MP1 = M + 1
      IF (Alpha .EQ. (-ONE)) THEN
        DO i = MP1,N,5
          X(i)     = -X(i)
          X(i + 1) = -X(i + 1)
          X(i + 2) = -X(i + 2)
          X(i + 3) = -X(i + 3)
          X(i + 4) = -X(i + 4)
        END DO
      ELSEIF (Alpha .EQ. ZERO) THEN
        DO i = MP1,N,5
          X(i)     = ZERO
          X(i + 1) = ZERO
          X(i + 2) = ZERO
          X(i + 3) = ZERO
          X(i + 4) = ZERO
        END DO
      ELSE
        DO i = MP1,N,5
          X(i)     = Alpha*X(i)
          X(i + 1) = Alpha*X(i + 1)
          X(i + 2) = Alpha*X(i + 2)
          X(i + 3) = Alpha*X(i + 3)
          X(i + 4) = Alpha*X(i + 4)
        END DO
      END IF

      END SUBROUTINE WSCAL

!--------------------------------------------------------------
      REAL(kind=dp) FUNCTION WLAMCH( C )
!--------------------------------------------------------------
!     returns epsilon machine
!     after LAPACK
!     replace this by the function from the optimized LAPACK implementation:
!          CALL SLAMCH('E') or CALL DLAMCH('E')
!--------------------------------------------------------------
!      USE aqchem_Precision

      CHARACTER ::  C
      INTEGER    :: i
      REAL(kind=dp), SAVE  ::  Eps
      REAL(kind=dp)  ::  Suma
      REAL(kind=dp), PARAMETER  ::  ONE=1.0_dp, HALF=0.5_dp
      LOGICAL, SAVE   ::  First=.TRUE.
      
      IF (First) THEN
        First = .FALSE.
        Eps = HALF**(16)
        DO i = 17, 80
          Eps = Eps*HALF
          CALL WLAMCH_ADD(ONE,Eps,Suma)
          IF (Suma.LE.ONE) GOTO 10
        END DO
        PRINT*,'ERROR IN WLAMCH. EPS < ',Eps
        RETURN
10      Eps = Eps*2
        i = i-1      
      END IF

      WLAMCH = Eps

      END FUNCTION WLAMCH
     
      SUBROUTINE WLAMCH_ADD( A, B, Suma )
!      USE aqchem_Precision
      
      REAL(kind=dp) A, B, Suma
      Suma = A + B

      END SUBROUTINE WLAMCH_ADD
!--------------------------------------------------------------


!--------------------------------------------------------------
      SUBROUTINE SET2ZERO(N,Y)
!--------------------------------------------------------------
!     copies zeros into the vector y:  y <- 0
!     after BLAS
!--------------------------------------------------------------
      
      INTEGER ::  i,M,MP1,N
      REAL(kind=dp) ::  Y(N)
      REAL(kind=dp), PARAMETER :: ZERO = 0.0d0

      IF (N.LE.0) RETURN

      M = MOD(N,8)
      IF( M .NE. 0 ) THEN
        DO i = 1,M
          Y(i) = ZERO
        END DO
        IF( N .LT. 8 ) RETURN
      END IF    
      MP1 = M+1
      DO i = MP1,N,8
        Y(i)     = ZERO
        Y(i + 1) = ZERO
        Y(i + 2) = ZERO
        Y(i + 3) = ZERO
        Y(i + 4) = ZERO
        Y(i + 5) = ZERO
        Y(i + 6) = ZERO
        Y(i + 7) = ZERO
      END DO

      END SUBROUTINE SET2ZERO


!--------------------------------------------------------------
      REAL(kind=dp) FUNCTION WDOT (N, DX, incX, DY, incY) 
!--------------------------------------------------------------
!     dot produce: wdot = x(1:N)*y(1:N) 
!     only for incX=incY=1
!     after BLAS
!     replace this by the function from the optimized BLAS implementation:
!         CALL SDOT(N,X,1,Y,1) or  CALL DDOT(N,X,1,Y,1)
!--------------------------------------------------------------
!      USE messy_mecca_kpp_Precision
!--------------------------------------------------------------
      IMPLICIT NONE
      INTEGER :: N, incX, incY
      REAL(kind=dp) :: DX(N), DY(N) 

      INTEGER :: i, IX, IY, M, MP1, NS
                                 
      WDOT = 0.0D0 
      IF (N .LE. 0) RETURN 
      IF (incX .EQ. incY) IF (incX-1) 5,20,60 
!                                                                       
!     Code for unequal or nonpositive increments.                       
!                                                                       
    5 IX = 1 
      IY = 1 
      IF (incX .LT. 0) IX = (-N+1)*incX + 1 
      IF (incY .LT. 0) IY = (-N+1)*incY + 1 
      DO i = 1,N 
        WDOT = WDOT + DX(IX)*DY(IY) 
        IX = IX + incX 
        IY = IY + incY 
      END DO 
      RETURN 
!                                                                       
!     Code for both increments equal to 1.                              
!                                                                       
!     Clean-up loop so remaining vector length is a multiple of 5.      
!                                                                       
   20 M = MOD(N,5) 
      IF (M .EQ. 0) GO TO 40 
      DO i = 1,M 
         WDOT = WDOT + DX(i)*DY(i) 
      END DO 
      IF (N .LT. 5) RETURN 
   40 MP1 = M + 1 
      DO i = MP1,N,5 
          WDOT = WDOT + DX(i)*DY(i) + DX(i+1)*DY(i+1) + DX(i+2)*DY(i+2) +  &
                   DX(i+3)*DY(i+3) + DX(i+4)*DY(i+4)                   
      END DO 
      RETURN 
!                                                                       
!     Code for equal, positive, non-unit increments.                    
!                                                                       
   60 NS = N*incX 
      DO i = 1,NS,incX 
        WDOT = WDOT + DX(i)*DY(i) 
      END DO 

      END FUNCTION WDOT                                          


!--------------------------------------------------------------
      SUBROUTINE WADD(N,X,Y,Z)
!--------------------------------------------------------------
!     adds two vectors: z <- x + y
!     BLAS - like
!--------------------------------------------------------------
!     USE aqchem_Precision
      
      INTEGER :: i, M, MP1, N
      REAL(kind=dp) :: X(N),Y(N),Z(N)

      IF (N.LE.0) RETURN

      M = MOD(N,5)
      IF( M /= 0 ) THEN
         DO i = 1,M
            Z(i) = X(i) + Y(i)
         END DO
         IF( N < 5 ) RETURN
      END IF    
      MP1 = M+1
      DO i = MP1,N,5
         Z(i)     = X(i)     + Y(i)
         Z(i + 1) = X(i + 1) + Y(i + 1)
         Z(i + 2) = X(i + 2) + Y(i + 2)
         Z(i + 3) = X(i + 3) + Y(i + 3)
         Z(i + 4) = X(i + 4) + Y(i + 4)
      END DO

      END SUBROUTINE WADD
      
      
      
!--------------------------------------------------------------
      SUBROUTINE WGEFA(N,A,Ipvt,info)
!--------------------------------------------------------------
!     WGEFA FACTORS THE MATRIX A (N,N) BY
!           GAUSS ELIMINATION WITH PARTIAL PIVOTING
!     LINPACK - LIKE 
!--------------------------------------------------------------
!
      INTEGER       :: N,Ipvt(N),info
      REAL(kind=dp) :: A(N,N)
      REAL(kind=dp) :: t, dmax, da
      INTEGER       :: j,k,l
      REAL(kind=dp), PARAMETER :: ZERO = 0.0, ONE = 1.0

      info = 0

size: IF (n > 1) THEN
      
col:  DO k = 1, n-1

!        find l = pivot index
!        l = idamax(n-k+1,A(k,k),1) + k - 1
         l = k; dmax = abs(A(k,k))
         DO j = k+1,n
            da = ABS(A(j,k))
            IF (da > dmax) THEN
              l = j; dmax = da
            END IF
         END DO
         Ipvt(k) = l

!        zero pivot implies this column already triangularized
         IF (ABS(A(l,k)) < TINY(ZERO)) THEN
            info = k
            return
         ELSE   
            IF (l /= k) THEN
               t = A(l,k); A(l,k) = A(k,k); A(k,k) = t
            END IF
            t = -ONE/A(k,k)
            CALL WSCAL(n-k,t,A(k+1,k),1)
            DO j = k+1, n
               t = A(l,j)
               IF (l /= k) THEN
                  A(l,j) = A(k,j); A(k,j) = t
               END IF
               CALL WAXPY(n-k,t,A(k+1,k),1,A(k+1,j),1)
            END DO         
         END IF
         
       END DO col
       
      END IF size
      
      Ipvt(N) = N
      IF (ABS(A(N,N)) == ZERO) info = N
      
      END SUBROUTINE WGEFA


!--------------------------------------------------------------
      SUBROUTINE WGESL(Trans,N,A,Ipvt,b)
!--------------------------------------------------------------
!     WGESL solves the system
!     a * x = b  or  trans(a) * x = b
!     using the factors computed by WGEFA.
!
!     Trans      = 'N'   to solve  A*x = b ,
!                = 'T'   to solve  transpose(A)*x = b
!     LINPACK - LIKE 
!--------------------------------------------------------------

      INTEGER       :: N,Ipvt(N)
      CHARACTER     :: trans
      REAL(kind=dp) :: A(N,N),b(N)
      REAL(kind=dp) :: t
      INTEGER       :: k,kb,l

      
      SELECT CASE (Trans)

      CASE ('n','N')  !  Solve  A * x = b

!        first solve  L*y = b
         IF (n >= 2) THEN
          DO k = 1, n-1
            l = Ipvt(k)
            t = b(l)
            IF (l /= k) THEN
               b(l) = b(k)
               b(k) = t
            END IF
            CALL WAXPY(n-k,t,a(k+1,k),1,b(k+1),1)
          END DO
         END IF
!        now solve  U*x = y
         DO kb = 1, n
            k = n + 1 - kb
            b(k) = b(k)/a(k,k)
            t = -b(k)
            CALL WAXPY(k-1,t,a(1,k),1,b(1),1)
         END DO
      
      CASE ('t','T')  !  Solve transpose(A) * x = b

!        first solve  trans(U)*y = b
         DO k = 1, n
            t = WDOT(k-1,a(1,k),1,b(1),1)
            b(k) = (b(k) - t)/a(k,k)
         END DO
!        now solve trans(L)*x = y
         IF (n >= 2) THEN
         DO kb = 1, n-1
            k = n - kb
            b(k) = b(k) + WDOT(n-k,a(k+1,k),1,b(k+1),1)
            l = Ipvt(k)
            IF (l /= k) THEN
               t = b(l); b(l) = b(k); b(k) = t
            END IF
         END DO
         END IF
   
      END SELECT

      END SUBROUTINE WGESL
! End of BLAS_UTIL function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



END MODULE aqchem_LinearAlgebra

